Mitochondrial complex 1 gene analysis in keratoconus.

PURPOSE
Keratoconus is characterized by the thinning of corneal stroma, resulting in reduced vision. The exact etiology of keratoconus (KC) is still unknown. The involvement of oxidative stress (OS) in this disease has been reported. However, the exact mechanism of OS in keratoconus is still unknown. Thus we planned this study to screen mitochondrial complex I genes for sequence changes in keratoconus patients and controls, as mitochondrial complex I is the chief source of reactive oxygen species (ROS) production.


METHODS
A total of 20 keratoconus cases and 20 healthy controls without any ocular disorder were enrolled in this study. Mitochondrial complex I genes (ND1, 2, 3, 4, 4L, 5, and 6) were amplified in all patients and controls using 12 pairs of primers by PCR. After sequencing, DNA sequences were analyzed against the mitochondrial reference sequence NC_012920. Haplogroup frequency based Principle Component Analysis (PCA) was constructed to determine whether the gene pool of keratoconus patients is closer to major populations in India.


RESULTS
DNA sequencing revealed a total 84 nucleotide variations in patients and 29 in controls. Of 84 nucleotide changes, 18 variations were non-synonymous and two novel frame-shift mutations were detected in cases. Non-synonymous mtDNA sequence variations may account for increased ROS and decreased ATP production. This ultimately leads to OS; which is a known cause for variety of corneal abnormalities. Haplotype analysis showed that most of the patients were clustered under the haplogroups: T, C4a2a, R2'TJ, M21'Q1a, M12'G2a2a, M8'CZ and M7a2a, which are present as negligible frequency in normal Indian population, whereas only few patients were found to be a part of the other haplogroups like U7 (Indo-European), R2 and R31, whose origin is contentious.


CONCLUSIONS
Mt complex I sequence variations are the main cause of elevated ROS production which leads oxidative stress. This oxidative stress then starts a cascade of events which ultimately can lead to keratoconus. Prompt antioxidant therapy should be initiated in keratoconus patients to minimize ROS related damage.

these patients, 14 were males and 6 were females. The mean age of presentation was 17.2 years. Diagnosis of keratoconus involved the presence of characteristic topographic features, such as inferior or central corneal steepening, or an asymmetric bowtie pattern with skewing of the radial axes, and the presence of one or more of the following characteristic, clinical features in one or both eyes: conical corneal deformation, munsen sign, corneal stromal thinning, a Fleischer ring or Vogt striae. Family history up to three generations was collected and pedigrees were drawn. All 20 cases were sporadic without any family history. All keratoconus cases secondary to causes like trauma, surgery, Ehlers Danlos syndrome, osteogenesis imperfecta, and pellucid marginal degeneration were excluded from the study. Twenty ethnically matched normal individuals without any ocular disorder were enrolled as controls. Health information was obtained from controls through the questionnaire; all underwent ophthalmological examination. Five milliliters of blood was collected by venipuncture in EDTA (EDTA) vaccutainers (Greiner Bio-One GmbH, Frickenhausen, Germany) from both patients and controls. DNA was extracted from whole blood samples by the inorganic method. For the population study, controls were taken from published data defining the lineage of the Indian population [14,15].

Polymerase chain reaction (PCR) amplification and sequence analysis of the mitochondrial DNA coding region:
The mitochondrial complex 1 (ND1, ND2, ND3, ND4, ND4L, ND5, and ND6 [ND stands for NADH dehydrogenase]) was amplified in all patients and controls using 12 pairs of primers using cycling conditions as described by Kumar and associates [16] and presented in Table 2. Briefly, PCR amplifications for all primer sets were performed in a 40-μl volume containing 1.0 μl of 20 μM stock solution for each primer (Eurofins Genomics India pvt Ltd, Bangalore, India), 100 ng of genomic DNA, 1 unit of Taq polymerase (Banglore Genei, Bengaluru, Karnataka, India), 0.1 mM of each deoxynucleotide triphosphate (dNTP), and 4 μl of 10× PCR buffer (with 15 mM MgCl2) by means of 30 cycles of amplification, each consisting of 30 s denaturation at 94 °C, 30 s annealing at 55 °C, and 1 min extension at 72 °C. Finally, an extension for 5 min at 72 °C was performed. Amplified PCR products were purified using a gel/PCR DNA fragments extraction kit (catalog number DF100; Geneaid Biotech Ltd., All sequence variants from both KC patients and controls were compared to the Human Mitochondrial reference sequence NC_012920 provided by the National Center for Biotechnology Information (NCBI) using ClustalW2 (multiple sequence alignment program for DNA; European Molecular Biology Laboratory (EMBL)-European Bioinformatics Institute (EBI).
Computational assessment of missense mutations: For prediction of pathogenic characteristics of all nonsynonymous mtDNA changes two homology based programs PolyPhen-2 (Polymorphism Phenotyping) and SIFT (Sorting PolyPhen structurally analyzes an amino acid polymorphism and predicts whether that amino acid change is likely to be deleterious to protein function [17][18][19]. Polyphen-2 is more advanced version of the earlier version PolyPhen [20]. The prediction is based on the positionspecific independent counts (PSIC) score derived from multiple sequence alignments of observations in case of functional domain of protein and predicted hydrophobic and transmembrane (PHAT) matrix element difference in case of transmembrane region of protein. PolyPhen scores of above 0.85 indicate the polymorphism is probably damaging to protein function. Scores of above 0.15 are possibly damaging, and scores of less than 0.15 are classified as benign.
SIFT is a sequence homology-based tool that sorts intolerant from tolerant amino acid substitutions and predicts whether an amino acid substitution in a protein will have a phenotypic effect [21][22][23]. SIFT is based on the premise that protein evolution is correlated with protein function. Positions important for function should be conserved in an alignment of the protein family, whereas unimportant positions should appear diverse in an alignment. Positions with normalized probabilities less than 0.05 are predicted to be deleterious and, those greater than or equal to 0.05 are predicted to be tolerated. Haplogroup and phylogenetic analysis: To check the fidelity of our conclusion, the evolutionary information and the significance of mutations should be known. For haplogroups (Hg) analysis we have carefully chosen two hundred healthy individual samples from same area for comparison analysis and these were also treated as controls. For all control samples, sequences of the control region were determined from position 16024 to 00300, using the ABI Prism Dye Terminator cycle-sequencing protocols developed by Applied Biosystems (Perkin-Elmer, Foster City, CA), to provide an initial presumed Hg assignment and cases were haplogrouped by complete coding region sequences. The C-track length variation at positions 16182 and 16183 in HVS-I and the indels at positions 00309 and 00315 in HVS-II were excluded from further analyses. Hg assignment was then confirmed, based on control and coding region Hg defining polymorphisms determined by means of direct sequencing.
The NETWORK 4.5.1.6 program was run for placing all the mutations of control samples in their respective phylogenetic tree using the protocol as described at the Fluxus Engineering Website. The matrilineal lineages of the case were drawn manually in the reduced median network of control samples, to create the topology map we have applied the reduced median algorithm (r=1), followed by the medianjoining algorithm (epsilon=2).  The MVSP software package (Kovach WL, Services KC. MVSP -A multi-variate statistical package for Windows ver 3.13m. 2004) was used to identify the principal components (PCs) of mitochondrial variations that lead to form a haplogroup for every individual. To express the relative importance of top two eigenvectors in the resulting PCA plot, two axes were scaled by a factor equal to the square root of the corresponding eigen value. This experiment was repeated to confirm the outcomes.

RESULTS
Sequence variation in Complex I genes: DNA sequencing of Complex I genes revealed a total 84 nucleotide variation in patients (Table 3) and 29 variations in controls (Table 4). Of the 84 nucleotide variations in patients, 18 (21.42%) were non synonymous (Table 5), 52 (61.90%) were synonymous, 9 (10.71%) variations were in RNA genes and 3 (3.57%) were Of 84 variations, 2 variations were frame-shift (11273G>A, 5300T>T). In one patient (KC 16) a single base deletion of guanine was observed at mtDNA position 11273. This caused a frame shift mutation after codon 172 (Gly>Ala) and introduced a stop codon at position 174 which resulted in a 173 amino acids truncated protein.This variation was homoplasmic ( Figure 1).
In patient KC 2, we found a 2 base pair (CA) deletion at genomic position 5300 and 5301. This frame shift mutation altered the amino acid reading frame in ND2 protein at position 277. This CA deletion produced a truncated protein of 287 amino acids ( Figure 2).
In silico analysis: SIFT analysis revealed two pathogenic changes (p.L71I and p.T9A) and PolyPhen revealed two pathogenic changes (p.N30L and p.I172V). The polyphen score of p.T9A was not available (no result for this mutation was available through PolyPhen; Table 5).  Principle component analysis: The tight cluster in Principal component Analysis (PCA) plot comprises the north-western, western, and north Indian population whereas the southern Indian and eastern Indian population is caught in a loose cluster (Figure 3). The controls were taken from published data defining Indian lineage for PCA and Haplogroup Network. We have treated the patients as a sub group of individuals having genetic structure different from normal Indians e.g., population. The inferences from PCA plot strongly supports our motive behind the planning of experiment, interestingly the patient population has not shown any relevant genetic affinity with other macropopulations of India.

DISCUSSION
In this study we analyzed mitochondrial complex 1 gene in 20 keratoconus patients (negative for VSX1 mutations [13]) and 20 unrelated healthy controls. The cornea, being an avascular structure and the first in line of ultraviolet (UV) radiation, is very susceptible to UV induced oxidative damage. Previous studies [10,11] suggested the role of oxidative stress in corneal disorders and congenital glaucoma [12]. Since the complex 1 NADH group of genes are most frequently associated with increased ROS production and oxidative stress [12,24], in this pilot study we analyzed mitochondrial complex I gene for sequence variations. Most of the mutations were found in ND5 (n=28) followed by ND4 (15) and then ND2 (13). The frequent variations in ND5 are in accordance with previous reports that mutations in ND5 gene of complex 1 play an important role in mitochondrial diseases [25].
In this study we report two novel frame shift mutations. Patient (KC 2) harbored a two base deletion (CA) which caused a frameshift and introduced a stop codon at position 287 in protein (normal ND2 protein is 347 amino acids long). The truncated protein cannot substitute the wild type ND2 protein as frameshift altered the reading frame of ND2. Sequence variations in this gene are associated with several diseases e.g., Leigh syndrome, breast cancer, myocardial infarction, Parkinson disease, and primary congenital glaucoma (PCG) [12,[26][27][28][29].
Patient (KC 16) harbored a single base deletion which resulted in a frame shift mutation after codon 172 (Gly>Ala) and introduced a stop codon at position 174 in protein and produced a truncated protein of 173 amino acids (wild type ND4 protein is 459 amino acids long).
Studies have documented that G10398A is associated with elevated ROS production due to altered complex 1 function [29][30][31][32]. Role of this allele G10398A has been implicated in diseases like congenital glaucoma, Parkinson, Type-2 diabetes, and in pre-term births [12,[29][30][31][32][33]. The G10398A variation though associated with high ROS levels was present significantly higher in cases as compared to controls however this is present in 43% Indian population. The 4216T>C variation considered as secondary or intermediate LHON-Leber's Hereditary optic neuropathy mutation was also present in 3 patients. However these patients had no features of LHON.

Evolutionary insight of Mt complex I sequence variations:
The genetic diversity in India is very complex. Several mutations from even control regions have been classified into the associative agent for various diseases [34]. The degree of haplotype sharing between populations is to investigate the combined frequency of the shared haplotypes in two population groups. Thus, among the northern and the southern population groups the combined frequency of the haplotypes present also in the other group is significantly lower than that which we observed in the case of random groups. This is not surprising because West Eurasian-specific mtDNA haplogroups are rather frequent in northwest India [35]. Because the Indo-European and the Dravidic speakers of India are largely concentrated to the northern and southern parts of the subcontinent, respectively, the differences arising from geographic division of the Indian populations also  correspond to different linguistic groupings [36]. In this study, we found that all the mutations were apparently North-Indian specific with some novel mutations. The sequencing of Complex 1 revealed 84 mutations, of which 14, including 2 frame shift mutations and 4 non-synonymous mutations, were novel and exclusively observed in KC patients. Interestingly, most of the patients and their maternal relatives were clustered under the haplogroups (T, C4a2a, R2'TJ, M21'Q1a, M12'G2a2a, M8'CZ, M7a2a, U5b1, U1a3) which are present as negligible frequency in normal Indian population, whereas only few patients were found to be a part of the haplogroups whose origin is contentious i.e., U7 (Indo-European), R2 and R31. We have found three patients who fall under Indian haplogroups (M4, M4'63, R31a1) but they also carry the same sets of novel synonymous and non-synonymous mutations i.e., 4769, 4985, 5580, and 12850 ( Figure 4). We have found some novel mutations in addition to each individual's lineages and they are different from each other. This finding suggests the positive/causative role of different combinations of the mitochondrial coding mutations in this disease, as the normal population, completely lack these mutations. The patients harbored some novel mutations at the different sites in mitochondria i.e., 4769, 4985, 5580, and 12850. These variants have never been reported in any of the population studies whereas they were present in every patient. Nevertheless, it is impossible from the evolutionary point of view that these sets of mutations in the individuals from different haplogroup. By keeping in mind about mutation rates in the coding region and its natural selection [37], we propose that these variants could theoretically influence the patient's phenotype. However, the variants present in coding regions of mitochondrial gene are not conserved in course of evolution. The patients were apparently homoplasmic (only one type of mtDNA was present). To determine whether the maternal inherited gene pool of keratoconus patients is truly closer to any major populations in India, we have constructed the haplogroup frequency based PCA plot for mtDNA ( Figure  3). Indeed, this analysis shows ambiguously that the three Indian populations clusters tightly among themselves viz. North, North West and West populations and two populations are to be found in a loose cluster viz South Indian and East Indian, whereas the keratoconus population matches with none of them in the mtDNA PCA plots. However the genetic data indicates that the keratoconus patients comprise several different haplotypes, if they are compared to normal populations around them. Most of the patients are in the clades which are nonspecific to Indian lineages. This information suggests that keratoconus patients are among those who are recent migrants into India and some genes in mitochondria have acquired mutations which are not filtered by purifying selection. Our results explain that the patients are genetically unrelated to each other due to the present maternal lineages which were diversified in the history of evolution. This fact suggests that the polymorphisms which are playing pivotal roles in causing the disease are recently accumulated in the mitochondrial coding regions of an individual patient. We have found that mutations specifically found in KC patients can affect transcription, translation or have synergistic effect with other variants in causing the disease. It has been reported many times about synergistic effect of different mutations in mitochondria that can cause many severe diseases [38]. Nevertheless, it is impossible from the evolutionary point of view these sets of mutations to occur in the normal individual from different haplogroups. By keeping in mind about mutation rates in the coding region and its natural selection, we propose that these variants could theoretically influence the patient's disease. Non-synonymous mutations and frame shift mutations adversely affect C1 synergetics resulting in increase ROS production and mitochondrial dysfunction. KC corneas are unable to process ROS due to depleted or low ATP levels and increased ROS production and thereby undergo oxidative damage. These corneas have increased levels of malondialdehyde (MDA), which can results in altered protein function leading to cascade of events, including apoptosis that can damage the corneal tissues.
Thus this pilot study highlights the role of sequence variation in mitochondrial complex I gene in keratoconus patients. Such cases with elevated free radicals levels and oxidative damage to cornea may benefit immensely by antioxidant therapy.